# Identity Motivated Reasoning

source("OtherScripts/IMR-2-4-19-recode.r")
options(contrasts=c("contr.treatment", "contr.treatment"))
library(mediation)
source("OtherScripts/IMR-2-4-19-alternatevpcoding.r")


vpaltlong <- data.frame(timing=as.vector(unlist(timing)), openlength=as.vector(unlist(openlength)), qnumber=paste("Q", as.vector(unlist(sapply(1:(dim(timing)[2]), function(x) rep(x, dim(timing)[1])))), sep=""), black=rep(vpalt$black, dim(timing)[2]), resp=paste("Resp", rep(1:length(potmedian), dim(timing)[2]), sep="_"))

vpaltlong$lt <- log(vpaltlong$timing)

vpaltlong$proofficer <- vpaltlong$qnumber %in% c("Q1","Q3","Q4","Q7")

library(lme4)

vpaltlforpred <- data.frame(black=c(FALSE, TRUE, FALSE, TRUE), proofficer=c(TRUE, TRUE, FALSE, FALSE), resp=NA, qnumber=NA)

vpaltlreg <- lmer(lt~black*proofficer+(1|resp)+(1|qnumber), data=vpaltlong)
vpaltopenreg <- lmer(openlength~black*proofficer+(1|resp)+(1|qnumber), data=vpaltlong)


predict.fun <- function(my.lmm) {
  predict(my.lmm, newdata = vpaltlforpred, re.form=NA)   # This is predict.merMod 
}

time.lmm.boots <- confint(bootMer(vpaltlreg, predict.fun, nsim = 100))
predict(vpaltlreg, newdata=vpaltlforpred, re.form=NA)

write.csv(summary(vpaltlreg)$coefficients, "Tables/NoCutoff/PredictingTimeCoefs.csv")

amt.lmm.boots <- confint(bootMer(vpaltopenreg, predict.fun, nsim = 100))
predict(vpaltopenreg, newdata=vpaltlforpred, re.form=NA)

write.csv(summary(vpaltopenreg)$coefficients, "Tables/NoCutoff/PredictingWordsCoefs.csv")

pdf("Figures/NoCutoff/MultilevelModelWordsTime.pdf", width=16, height=8)
par(mfrow=c(1,2))
plot(c(.5, 2.5), c(log(20), log(200)), type="n", axes=FALSE, xlab="Statement Direction", ylab="Time Spent Reading (Seconds, Log Scale)", main="a) Predicting Time Spent Reading")
axis(1, c(1,2), labels=c("Pro-Victim", "Pro-Officer"))
axis(1, c(-99,99))
axis(2, log(seq(20, 400, 20)), seq(20, 400, 20))
abline(h=log(seq(20, 400, 20)), lty=3, col="gray")
axis(2, c(-99,9999))
axis(3, c(-99,9999))
axis(4, c(-99,9999))
segments(c(1.9,2.1,.9,1.1), time.lmm.boots[,1], c(1.9,2.1,.9,1.1), time.lmm.boots[,2], col=c("dark gray", "black", "dark gray", "black"), lwd=2)
lines(c(1.9,2.1,.9,1.1), predict(vpaltlreg, newdata=vpaltlforpred, re.form=NA), col=c("dark gray", "black", "dark gray", "black"), pch=20, type="p", lwd=3)
plot(c(.5, 2.5), c(20, 140), type="n", axes=FALSE, xlab="Statement Direction", ylab="Characters Written", main="b) Predicting Characters Written")
axis(1, c(1,2), labels=c("Pro-Victim", "Pro-Officer"))
axis(1, c(-99,99))
axis(2, seq(20, 400, 20), seq(20, 400, 20))
abline(h=seq(20, 400, 20), lty=3, col="gray")
axis(2, c(-99,9999))
axis(3, c(-99,9999))
axis(4, c(-99,9999))
segments(c(1.9,2.1,.9,1.1), amt.lmm.boots[,1], c(1.9,2.1,.9,1.1), amt.lmm.boots[,2], col=c("dark gray", "black", "dark gray", "black"), lwd=2)
lines(c(1.9,2.1,.9,1.1), predict(vpaltopenreg, newdata=vpaltlforpred, re.form=NA), col=c("dark gray", "black", "dark gray", "black"), pch=20, type="p", lwd=3)
dev.off()

# Table 2 - Perceptions by Respondent Race

pnames <- c("Officer's Actions Appropriate", "Officer Should Be Charged", "Likely That Victim Attacked Officer", "Victim Had Weapon", "Role of Race In Shooting")

#p01 <- as.data.frame(sapply(principals, nalevs))

pdifs <- with(vpalt, t(sapply(p01, function(x) unlist(wtd.t.test(nalevs(x)[black==FALSE], nalevs(x)[black==TRUE])))))

whitemnsd <- with(vpalt, t(sapply(p01, function(x) c(M=mean(nalevs(x)[black==FALSE], na.rm=TRUE), SD=sd(nalevs(x)[black==FALSE], na.rm=TRUE)/sqrt(length(!is.na(nalevs(x)[black==FALSE])))))))
blackmnsd <- with(vpalt, t(sapply(p01, function(x) c(M=mean(nalevs(x)[black==TRUE], na.rm=TRUE), SD=sd(nalevs(x)[black==TRUE], na.rm=TRUE)/sqrt(length(!is.na(nalevs(x)[black==TRUE])))))))

blkmin <- blackmnsd[,1]-1.96*blackmnsd[,2]
blkmax <- blackmnsd[,1]+1.96*blackmnsd[,2]
whtmin <- whitemnsd[,1]-1.96*whitemnsd[,2]
whtmax <- whitemnsd[,1]+1.96*whitemnsd[,2]

pdifsout <- as.data.frame(pdifs, stringsAsFactors=FALSE)
pdifsout$coefficients.p.value.adjusted <- p.adjust(as.numeric(as.character(pdifsout$coefficients.p.value)), method="BY")

pdifstable <- data.frame(Whites=pdifs[,"additional.Mean.x"], Blacks=pdifs[,"additional.Mean.y"], Difference=pdifs[,"additional.Difference"], sig=starmaker(as.numeric(as.character(pdifs[,"coefficients.p.value"]))))

rownames(pdifstable) <- pnames

pdifstableNs <-  t(sapply(p01, function(x) c(White=sum(!is.na(x) & vpalt$black==FALSE), Black=sum(!is.na(x) & vpalt$black==TRUE))))

write.csv(pdifstable, "Tables/TableH2-PerceptionsOfNewIncidentByRacialCategory.csv")

pdifsappendix <- data.frame(pdifsout[,c("additional.Mean.x", "additional.Mean.y", "additional.Difference", "additional.Std. Err", "coefficients.t.value", "coefficients.p.value", "coefficients.p.value.adjusted")], pdifstableNs)
colnames(pdifsappendix) <- c("Whites Mean", "Blacks Mean", "Difference", "Standard Error", "t-statistic", "p-value", "adjusted p-value", "N of Whites", "N of Blacks")
rownames(pdifsappendix) <- pnames

write.csv(pdifsappendix, "Tables/TableHG1-FullInformationOnPerceptionsOfNewIncidentByRacialCategory.csv")

# Regressions Predicting Perceptions

ext <- function(x)
    cbind(Coef=coef(summary(x))[,"Estimate"], SE=coef(summary(x))[,"Std. Error"], p=coef(summary(x))[,(colnames(coef(summary(x))) %in% c("Pr(>|t|)", "Pr(>|z|)"))], padj=p.adjust(coef(summary(x))[,(colnames(coef(summary(x))) %in% c("Pr(>|t|)", "Pr(>|z|)"))], "BY"))

p01regs <- with(vpalt, lapply(p01, function(y) lm(y~black+female+age01+ed01+pid01+libcon01)))

p01regsNs <- sapply(p01regs, function(x) length(predict(x)))
p01regsR2s <- sapply(p01regs, function(x) summary(x)$r.squared)

p01regsout <- rbind(pnames[1], ext(p01regs[[1]]), N=p01regsNs[1], R2=p01regsR2s[1])
for(i in 2:length(p01regs)){
    regtemp <- rbind(pnames[i], ext(p01regs[[i]]), N=p01regsNs[i], R2=p01regsR2s[i])
    p01regsout <- cbind(p01regsout, regtemp)
}

matrix(p.adjust(as.vector(p01regsout[!(rownames(p01regsout) %in% c("", "N", "R2")), colnames(p01regsout)=="p"]), "BY"), ncol=i, byrow=FALSE)

p01regsout[!(rownames(p01regsout) %in% c("", "N", "R2")), colnames(p01regsout)=="padj"]

write.csv(p01regsout, "Tables/TableHG2-RegressionsPredictingPerceptions.csv")

## Statement Evaluations by Respondent Race

snames <- c("Police Chief (Pro-Officer)", "Mr. Davis (Pro-Victim)", "Officer Silver (Pro-Officer)", "Mrs. Walker (Pro-Officer)", "Mrs. Thomas (Pro-Victim)", "Mrs. Williams (Pro-Victim)", "Mr. Anthony (Pro-Officer)")

snamesplus <- c(snames, "All Pro-Officer", "All Pro-Victim", "Pro Officer -Walker")

w01 <- as.data.frame(sapply(weightset, nalevs))
a01 <- as.data.frame(sapply(accurateset, nalevs))
b01 <- as.data.frame(sapply(biasedset, nalevs))
ec01 <- as.data.frame((w01+a01+(1-b01)))/3

w01$allprooff <- rowMeans(w01[c(1,3,4,7)])
w01$allprotaylor <- rowMeans(w01[c(2,5,6)])
w01$allprooffNoWalker <- rowMeans(w01[c(1,3,7)])

a01$allprooff <- rowMeans(a01[c(1,3,4,7)])
a01$allprotaylor <- rowMeans(a01[c(2,5,6)])
a01$allprooffNoWalker <- rowMeans(a01[c(1,3,7)])

b01$allprooff <- rowMeans(b01[c(1,3,4,7)])
b01$allprotaylor <- rowMeans(b01[c(2,5,6)])
b01$allprooffNoWalker <- rowMeans(b01[c(1,3,7)])

ec01$allprooff <- rowMeans(ec01[c(1,3,4,7)])
ec01$allprotaylor <- rowMeans(ec01[c(2,5,6)])
ec01$allprooffNoWalker <- rowMeans(ec01[c(1,3,7)])

wdifs <- with(vpalt, t(sapply(w01, function(x) unlist(wtd.t.test(nalevs(x)[black==FALSE], nalevs(x)[black==TRUE])))))
adifs <- with(vpalt, t(sapply(a01, function(x) unlist(wtd.t.test(nalevs(x)[black==FALSE], nalevs(x)[black==TRUE])))))
bdifs <- with(vpalt, t(sapply(b01, function(x) unlist(wtd.t.test(nalevs(x)[black==FALSE], nalevs(x)[black==TRUE])))))
ecdifs <- with(vpalt, t(sapply(ec01, function(x) unlist(wtd.t.test(nalevs(x)[black==FALSE], nalevs(x)[black==TRUE])))))

wdifspadj <- p.adjust(wdifs[,"coefficients.p.value"], "BY")
adifspadj <- p.adjust(adifs[,"coefficients.p.value"], "BY")
bdifspadj <- p.adjust(bdifs[,"coefficients.p.value"], "BY")
ecdifspadj <- p.adjust(ecdifs[,"coefficients.p.value"], "BY")

fulltableadj <- p.adjust(c(wdifs[,"coefficients.p.value"], adifs[,"coefficients.p.value"], bdifs[,"coefficients.p.value"]), "BY")#, ecdifs[,"coefficients.p.value"]

wdifstable <- data.frame(Whites=wdifs[,"additional.Mean.x"], Blacks=wdifs[,"additional.Mean.y"], Difference=wdifs[,"additional.Difference"], sig=starmaker(as.numeric(as.character(wdifs[,"coefficients.p.value"]))))
adifstable <- data.frame(Whites=adifs[,"additional.Mean.x"], Blacks=adifs[,"additional.Mean.y"], Difference=adifs[,"additional.Difference"], sig=starmaker(as.numeric(as.character(adifs[,"coefficients.p.value"]))))
bdifstable <- data.frame(Whites=bdifs[,"additional.Mean.x"], Blacks=bdifs[,"additional.Mean.y"], Difference=bdifs[,"additional.Difference"], sig=starmaker(as.numeric(as.character(bdifs[,"coefficients.p.value"]))))
ecdifstable <- data.frame(Whites=ecdifs[,"additional.Mean.x"], Blacks=ecdifs[,"additional.Mean.y"], Difference=ecdifs[,"additional.Difference"], sig=starmaker(as.numeric(as.character(ecdifs[,"coefficients.p.value"]))))

rownames(wdifstable) <- rownames(adifstable) <- rownames(bdifstable) <- rownames(ecdifstable) <- snamesplus

sdifstable <- data.frame(Weight=wdifstable, Accuracy=adifstable, Biased=bdifstable, Combined=ecdifstable)
sdifs <- cbind(rbind(cbind(wdifs, wdifspadj), cbind(adifs, adifspadj), cbind(bdifs, bdifspadj)), fulltableadj)#, cbind(ecdifs, ecdifspadj)

write.csv(sdifstable[,1:4], "Tables/TableH3-WeightsByRace.csv")

Nws <-  sapply(w01, function(x) c(White=sum(!is.na(x) & vpalt$black==FALSE), Black=sum(!is.na(x) & vpalt$black==TRUE)))
Nas <-  sapply(a01, function(x) c(White=sum(!is.na(x) & vpalt$black==FALSE), Black=sum(!is.na(x) & vpalt$black==TRUE)))
Nbs <-  sapply(b01, function(x) c(White=sum(!is.na(x) & vpalt$black==FALSE), Black=sum(!is.na(x) & vpalt$black==TRUE)))


sdifsappendix <- data.frame(sdifs[,c("additional.Mean.x", "additional.Mean.y", "additional.Difference", "additional.Std. Err", "coefficients.t.value", "coefficients.p.value", "wdifspadj")], rbind(t(Nws), t(Nas), t(Nbs)))
colnames(sdifsappendix) <- c("Whites Mean", "Blacks Mean", "Difference", "Standard Error", "t-statistic", "p-value", "adjusted p-value (per outcome type)", "N of Whites", "N of Blacks")
rownames(sdifsappendix) <- as.vector(unlist(sapply(c("Weight", "Accuracy", "Bias"), function(x) paste(snamesplus, x))))

write.csv(sdifsappendix, "Tables/TableHG3-WeightsAccuracyBiasByRace.csv")


# Regressions Predicting Weights

w01regs <- with(vpalt, lapply(w01, function(y) lm(y~black+female+age01+ed01+pid01+libcon01)))

w01regsNs <- sapply(w01regs, function(x) length(predict(x)))
w01regsR2s <- sapply(w01regs, function(x) summary(x)$r.squared)

w01regsout <- rbind(snamesplus[1], ext(w01regs[[1]]), N=w01regsNs[1], R2=w01regsR2s[1])
for(i in 2:length(w01regs)){
    regtemp <- rbind(snamesplus[i], ext(w01regs[[i]]), N=w01regsNs[i], R2=w01regsR2s[i])
    w01regsout <- cbind(w01regsout, regtemp)
}

write.csv(w01regsout, "Tables/TableHG4-RegressionsPredictingWeights.csv")


# Correlations Between Statement Weights and Perceptions

corset <- wtd.cor(p01, w01)
corsetpadj <- matrix(p.adjust(corset$p.value, "BY"), nrow=dim(corset$p.value)[1], ncol=dim(corset$p.value)[2])
corout <- sapply(1:dim(p01)[2], function(x) paste(rd(corset$correlation[,x]), starmaker(corset$p.value[,x]), sep=""))
rownames(corout) <- rownames(corsetpadj) <- snamesplus
colnames(corout) <- colnames(corsetpadj) <- pnames

#write.csv(corsetpadj, "Tables/TableH4adjustedpvalues.csv")

write.csv(corout, "Tables/TableHG5-CorrelationsBtwStatementWeightsand RespondentPerceptions.csv")



# Regressions Predicting Perceptions Controlling for Expectations and Racial Identification

fd <- with(vpalt, data.frame(p01plus, sexfac, ed01, age01, pid01, libcon01, black, racres, differentialperceptions, salience, libcon01))
usable <- apply(fd, 1, function(x) sum(is.na(x))==0)

racebaseregs <- lapply(p01plus[usable,], function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+pid01+libcon01+black)))
racesalienceregs <- lapply(p01plus[usable,], function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+pid01+libcon01+black+salience*black)))
racepriorregs <- lapply(p01plus[usable,], function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+pid01+libcon01+black+racres+differentialperceptions)))
racefullregs <- lapply(p01plus[usable,], function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+pid01+libcon01+black+racres+differentialperceptions+salience*black)))

nracebaseregs <- sapply(racebaseregs, function(x) length(predict(x)))
nracesalienceregs <- sapply(racesalienceregs, function(x) length(predict(x)))
nracepriorregs <- sapply(racepriorregs, function(x) length(predict(x)))
nracefullregs <- sapply(racefullregs, function(x) length(predict(x)))

R2racebaseregs <- sapply(racebaseregs, function(x) summary(x)$r.squared)
R2racesalienceregs <- sapply(racesalienceregs, function(x) summary(x)$r.squared)
R2racepriorregs <- sapply(racepriorregs, function(x) summary(x)$r.squared)
R2racefullregs <- sapply(racefullregs, function(x) summary(x)$r.squared)

salienceanova <- lapply(1:length(racebaseregs), function(x) anova(racebaseregs[[x]], racesalienceregs[[x]]))
priorsanova <- lapply(1:length(racebaseregs), function(x) anova(racebaseregs[[x]], racepriorregs[[x]]))
allvnoneanova <- lapply(1:length(racebaseregs), function(x) anova(racebaseregs[[x]], racefullregs[[x]]))
allvsalienceanova <- lapply(1:length(racebaseregs), function(x) anova(racesalienceregs[[x]], racefullregs[[x]]))
allvpaltriorsanova <- lapply(1:length(racebaseregs), function(x) anova(racepriorregs[[x]], racefullregs[[x]]))

R2table <- cbind(R2racebaseregs, R2racesalienceregs, SalienceImprovement=(R2racesalienceregs-R2racebaseregs), R2racepriorregs, PriorsImprovement=(R2racepriorregs-R2racebaseregs), R2racefullregs, PriorsDifference=(R2racefullregs-R2racesalienceregs), SalienceDifference=(R2racefullregs-R2racepriorregs))

write.csv(R2table, "Tables/TableHG6R2s.csv")

for(i in 1:length(racebaseregs)){
    coefcomb <- cbind(ext(racebaseregs[[i]])[match(names(coef(racefullregs[[i]])), names(coef(racebaseregs[[i]]))),], ext(racesalienceregs[[i]])[match(names(coef(racefullregs[[i]])), names(coef(racesalienceregs[[i]]))),], ext(racepriorregs[[i]])[match(names(coef(racefullregs[[i]])), names(coef(racepriorregs[[i]]))),], ext(racefullregs[[i]]))
    rownames(coefcomb) <- names(coef(racefullregs[[i]]))
    nvals <- c(nracebaseregs[i], "", "", "", nracesalienceregs[i], "", "", "", nracepriorregs[i], "", "", "", nracefullregs[i], "", "", "")
    R2vals <- c(R2racebaseregs[i], "", "", "", R2racesalienceregs[i], "", "", "", R2racepriorregs[i], "", "", "", R2racefullregs[i], "", "", "")
    outsetmod <- rbind(coefcomb, "", nvals, R2vals)
    write.csv(outsetmod, paste("Tables/TableHG6", LETTERS[i], "_RegressionsFor_", names(p01)[i], "_ControllingForExpectationsSalience.csv", sep=""))
}

anovasets <- NULL
for(i in 1:length(salienceanova))
    anovasets[[i]] <- as.data.frame(rbind(salienceanova[[i]][2,], priorsanova[[i]][2,], allvnoneanova[[i]][2,], allvsalienceanova[[i]][2,], allvpaltriorsanova[[i]][2,]))

anovaouts <- sapply(anovasets, function(x) paste("F = ", rd(x$F, 1), " df = ", x$Df, " p = ", rd(x[,'Pr(>F)']), sep=""))

rownames(anovaouts) <- c("Salience", "Priors", "Both vs None", "Both vs Salience", "Both vs Priors")
colnames(anovaouts) <- names(racebaseregs)

write.csv(anovaouts, "Tables/TableHG6Anovas.csv")


p01regsfull <- with(vpalt, lapply(p01, function(y) lm(y~black*salience+female+age01+ed01+pid01+libcon01)))

w01regsNs <- sapply(w01regs, function(x) length(predict(x)))
w01regsR2s <- sapply(w01regs, function(x) summary(x)$r.squared)

w01regsout <- rbind(snamesplus[1], ext(w01regs[[1]]), N=w01regsNs[1], R2=w01regsR2s[1])
for(i in 2:length(w01regs)){
    regtemp <- rbind(snamesplus[i], ext(w01regs[[i]]), N=w01regsNs[i], R2=w01regsR2s[i])
    w01regsout <- cbind(w01regsout, regtemp)
}

write.csv(w01regsout, "Tables/TableHG4-RegressionsPredictingWeights.csv")



# Summary Judgments When Identity is Primed or by Moderators

baseregs <- lapply(p01, function(y) with(vpalt, lm(y~condition+black+female+age01+ed01+pid01+libcon01)))
baseregsimp <- lapply(p01, function(y) with(vpalt, lm(y~mi01+black+female+age01+ed01+pid01+libcon01)))
baseregsclose <- lapply(p01, function(y) with(vpalt, lm(y~mc01+black+female+age01+ed01+pid01+libcon01)))
baseregslinked <- lapply(p01, function(y) with(vpalt, lm(y~ml01+black+female+age01+ed01+pid01+libcon01)))
baseregssalience <- lapply(p01, function(y) with(vpalt, lm(y~salience+black+female+age01+ed01+pid01+libcon01)))

conditionregs <- lapply(p01, function(y) with(vpalt, lm(y~condition*black+female+age01+ed01+pid01+libcon01)))
conditionregsIRplus <- lapply(p01, function(y) with(vpalt, lm(y~condition*salience*black+female+age01+ed01+pid01+libcon01)))
conditionregsIRminus <- lapply(p01, function(y) with(vpalt, lm(y~(condition+salience)*black+female+age01+ed01+pid01+libcon01)))

ARset <- lapply(1:length(conditionregsIRplus), function(x) anova(conditionregsIRplus[[x]], conditionregsIRminus[[x]]))


importanceregs <- lapply(p01, function(y) with(vpalt, lm(y~mi01*black+female+age01+ed01+pid01+libcon01)))
closenessregs <- lapply(p01, function(y) with(vpalt, lm(y~mc01*black+female+age01+ed01+pid01+libcon01)))
linkedregs <- lapply(p01, function(y) with(vpalt, lm(y~ml01*black+female+age01+ed01+pid01+libcon01)))
salienceregs <- lapply(p01, function(y) with(vpalt, lm(y~salience*black+female+age01+ed01+pid01+libcon01)))

ftests <- lapply(1:length(baseregs), function(x) anova(baseregs[[x]], conditionregs[[x]]))
fouts <- as.matrix(t(sapply(ftests, function(x) x[2,])))
rownames(fouts) <- pnames
fconds <- paste(rd(unlist(fouts[,"F"])), "(", unlist(fouts[,"Df"]), ",", unlist(fouts[,"Res.Df"]), ")", starmaker(unlist(fouts[,"Pr(>F)"])), sep="")
fcondsadj <- paste("F=", rd(unlist(fouts[,"F"])), "(", unlist(fouts[,"Df"]), ",", unlist(fouts[,"Res.Df"]), "); p=", rd(unlist(fouts[,"Pr(>F)"])), "; Adj.p=", rd(p.adjust(unlist(fouts[,"Pr(>F)"]), "BY")), sep="")

ftestsimp <- lapply(1:length(baseregs), function(x) anova(baseregsimp[[x]], importanceregs[[x]]))
foutsimp <- as.matrix(t(sapply(ftestsimp, function(x) x[2,])))
fcondsimp <- paste(rd(unlist(foutsimp[,"F"])), "(", unlist(foutsimp[,"Df"]), ",", unlist(foutsimp[,"Res.Df"]), ")", starmaker(unlist(foutsimp[,"Pr(>F)"])), sep="")
fcondsimpadj <- paste("F=", rd(unlist(foutsimp[,"F"])), "(", unlist(foutsimp[,"Df"]), ",", unlist(foutsimp[,"Res.Df"]), "); p=", rd(unlist(foutsimp[,"Pr(>F)"])), "; Adj.p=", rd(p.adjust(unlist(foutsimp[,"Pr(>F)"]), "BY")), sep="")

ftestsclose <- lapply(1:length(baseregs), function(x) anova(baseregsclose[[x]], closenessregs[[x]]))
foutsclose <- as.matrix(t(sapply(ftestsclose, function(x) x[2,])))
fcondsclose <- paste(rd(unlist(foutsclose[,"F"])), "(", unlist(foutsclose[,"Df"]), ",", unlist(foutsclose[,"Res.Df"]), ")", starmaker(unlist(foutsclose[,"Pr(>F)"])), sep="")
fcondscloseadj <- paste("F=", rd(unlist(foutsclose[,"F"])), "(", unlist(foutsclose[,"Df"]), ",", unlist(foutsclose[,"Res.Df"]), "); p=", rd(unlist(foutsclose[,"Pr(>F)"])), "; Adj.p=", rd(p.adjust(unlist(foutsclose[,"Pr(>F)"]), "BY")), sep="")

ftestslinked <- lapply(1:length(baseregs), function(x) anova(baseregslinked[[x]], linkedregs[[x]]))
foutslinked <- as.matrix(t(sapply(ftestslinked, function(x) x[2,])))
fcondslinked <- paste(rd(unlist(foutslinked[,"F"])), "(", unlist(foutslinked[,"Df"]), ",", unlist(foutslinked[,"Res.Df"]), ")", starmaker(unlist(foutslinked[,"Pr(>F)"])), sep="")
fcondslinkedadj <- paste("F=", rd(unlist(foutslinked[,"F"])), "(", unlist(foutslinked[,"Df"]), ",", unlist(foutslinked[,"Res.Df"]), "); p=", rd(unlist(foutslinked[,"Pr(>F)"])), "; Adj.p=", rd(p.adjust(unlist(foutslinked[,"Pr(>F)"]), "BY")), sep="")

ftestssalience <- lapply(1:length(baseregs), function(x) anova(baseregssalience[[x]], salienceregs[[x]]))
foutssalience <- as.matrix(t(sapply(ftestssalience, function(x) x[2,])))
fcondssalience <- paste(rd(unlist(foutssalience[,"F"])), "(", unlist(foutssalience[,"Df"]), ",", unlist(foutssalience[,"Res.Df"]), ")", starmaker(unlist(foutssalience[,"Pr(>F)"])), sep="")
fcondssalienceadj <- paste("F=", rd(unlist(foutssalience[,"F"])), "(", unlist(foutssalience[,"Df"]), ",", unlist(foutssalience[,"Res.Df"]), "); p=", rd(unlist(foutssalience[,"Pr(>F)"])), "; Adj.p=", rd(p.adjust(unlist(foutssalience[,"Pr(>F)"]), "BY")), sep="")

pvalsfouts <- unlist(fouts[,"Pr(>F)"])
pvalsfoutsimp <- unlist(foutsimp[,"Pr(>F)"])
pvalsfoutsclose <- unlist(foutsclose[,"Pr(>F)"])
pvalsfoutslinked <- unlist(foutslinked[,"Pr(>F)"])
pvalsfoutssalience <- unlist(foutssalience[,"Pr(>F)"])

pvalsfsset <- cbind(pvalsfouts, pvalsfoutsimp, pvalsfoutsclose, pvalsfoutslinked, pvalsfoutssalience)
pvalsfssetadj <- apply(pvalsfsset, 2, function(x) p.adjust(x, "BY"))
rownames(pvalsfssetadj) <- rownames(pvalsfsset)
colnames(pvalsfssetadj) <- colnames(pvalsfsset)

tab5 <- data.frame('Racial Primes'=fconds, 'ID Strength'=fcondsimp, 'Group Closeness'=fcondsclose, 'Linked Fate'=fcondslinked, 'All Salience'=fcondssalience)
tab5b <- data.frame('Racial Primes'=fcondsadj, 'ID Strength'=fcondsimpadj, 'Group Closeness'=fcondscloseadj, 'Linked Fate'=fcondslinkedadj, 'All Salience'=fcondssalienceadj)

rownames(tab5) <- pnames

#write.csv(tab5, "Tables/TableH5.csv")

#write.csv(pvalsfsset, "Tables/TableH5ps.csv")
#write.csv(pvalsfssetadj, "Tables/TableH5psadj.csv")

outcomeopts <- lapply(principals, function(x) gsub("No at", "No role at", gsub("NOT", "not", gsub("An e", "E", gsub("A m", "M", gsub("A l", "L", gsub(" have a weapon", "", gsub(" did", "", gsub(" should be charged", "", gsub("should NOT be charged", "not", gsub("He p", "P", gsub("He d", "D", gsub(" role", "", gsub(" likely", "", gsub(" appropriate", "", levels(x))))))))))))))))

pdf("Figures/FigureH1-ConditionandSalienceByOutcome.pdf", width=10, height=20)
par(mfrow=c(5,2))
par(oma=c(0,2,2,0))
for(i in 1:length(conditionregs)){
    plot(c(.5,3.5), 0:1, type="n", ylab="", ylim=0:1, xlab="Racial Priming Condition", main=paste(letters[2*(i-1)+1], ") Priming Condition\nF=", tab5[i,1], sep=""), axes=FALSE)
    axis(1, at=1:3, c("Control", "Identity", "ID + Exp"))
    axis(2, at=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), outcomeopts[[i]])
    axis(1, at=c(-99,99))
    axis(2, at=c(-99,99))
    axis(3, at=c(-99,99))
    axis(4, at=c(-99,99))
    abline(h=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), col="light gray", lty=3)
    #plotwtdinteraction(conditionregs[[i]], "condition", "black", bylevnames=c("White", "Black"), placement="topleft", add=TRUE, linecol=c("dark gray", "black"), secol="gray")
    fwc <- findwtdinteraction(conditionregs[[i]], "condition", "black", bylevnames=c("White", "Black"))#, placement="topleft", add=TRUE, linecol=c("dark gray", "black"), secol="gray")
    segments(1:length(fwc$Means$All["White",])-.1, fwc$Means$All["White",]+1.96*fwc$SEs$All["White",], 1:length(fwc$Means$All["White",])-.1, fwc$Means$All["White",]-1.96*fwc$SEs$All["White",], lwd=2, col="dark gray")
    segments(1:length(fwc$Means$All["Black",])+.1, fwc$Means$All["Black",]+1.96*fwc$SEs$All["Black",], 1:length(fwc$Means$All["Black",])+.1, fwc$Means$All["Black",]-1.96*fwc$SEs$All["Black",], lwd=2, col="black")
    lines(1:length(fwc$Means$All["White",])-.1, fwc$Means$All["White",], lwd=3, col="dark gray", pch=15, type="p")
    lines(1:length(fwc$Means$All["Black",])+.1, fwc$Means$All["Black",], lwd=3, col="black", pch=15, type="p")
    legend("topleft", legend=c("White", "Black"), col=c("dark gray", "black"), pch=15, lwd=2)
    mtext(pnames[i], 2, line=3)#, outer=TRUE)
    plot(0:1, 0:1, type="n", ylab="", ylim=0:1, xlab="Group Salience", main=paste(letters[2*(i-1)+2], ") Group Salience\nF=", tab5[i,5], sep=""), axes=FALSE)
    axis(1, at=seq(0, 1, 1/4), c("Lowest", "", "", "", "Highest"))
    axis(2, at=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), outcomeopts[[i]])
    axis(1, at=c(-99,99))
    axis(2, at=c(-99,99))
    axis(3, at=c(-99,99))
    axis(4, at=c(-99,99))
    abline(h=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), col="light gray", lty=3)
    plotwtdinteraction(salienceregs[[i]], "salience", "black", bylevnames=c("White", "Black"), placement="topleft", add=TRUE, linecol=c("dark gray", "black"), secol="gray")
}
mtext("Perceptions by Race and", 3, outer=TRUE)
dev.off()


pdf("Figures/FigureHG1-5by5_SeparateMeasures.pdf", width=14, height=20)
par(mfrow=c(5,5))
par(oma=c(0,2,2,0))
for(i in 1:length(conditionregs)){
    plot(c(1,3), 0:1, type="n", ylab="", ylim=0:1, xlab="Racial Priming Condition", main=paste(letters[5*(i-1)+1], ") Priming Condition\n", tab5b[i,1], sep=""), axes=FALSE)
    axis(1, at=1:3, c("Control", "Identity", "ID + Exp"))
    axis(2, at=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), outcomeopts[[i]])
    axis(1, at=c(-99,99))
    axis(2, at=c(-99,99))
    axis(3, at=c(-99,99))
    axis(4, at=c(-99,99))
    abline(h=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), col="light gray", lty=3)
    plotwtdinteraction(conditionregs[[i]], "condition", "black", bylevnames=c("White", "Black"), placement="topleft", add=TRUE, linecol=c("dark gray", "black"), secol="gray")
    mtext(pnames[i], 2, line=3)#, outer=TRUE)
    plot(0:1, 0:1, type="n", ylab="", ylim=0:1, xlab="Identity Importance", main=paste(letters[5*(i-1)+2], ") Identity Strength\n", tab5b[i,2], sep=""), axes=FALSE)
    axis(1, at=seq(0, 1, 1/4), c("Not at all", "A little", "Moderately", "Very", "Extremely"))
    axis(2, at=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), outcomeopts[[i]])
    axis(1, at=c(-99,99))
    axis(2, at=c(-99,99))
    axis(3, at=c(-99,99))
    axis(4, at=c(-99,99))
    abline(h=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), col="light gray", lty=3)
    plotwtdinteraction(importanceregs[[i]], "mi01", "black", bylevnames=c("White", "Black"), placement="topleft", add=TRUE, linecol=c("dark gray", "black"), secol="gray")
    plot(0:1, 0:1, type="n", ylab="", ylim=0:1, xlab="Closeness to Racial Group", main=paste(letters[5*(i-1)+3], ") Group Closeness\n", tab5b[i,3], sep=""), axes=FALSE)
    axis(1, at=seq(0, 1, 1/4), c("Not at all", "A little", "Somewhat", "Very", "Extremely"))
    axis(2, at=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), outcomeopts[[i]])
    axis(1, at=c(-99,99))
    axis(2, at=c(-99,99))
    axis(3, at=c(-99,99))
    axis(4, at=c(-99,99))
    abline(h=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), col="light gray", lty=3)
    plotwtdinteraction(closenessregs[[i]], "mc01", "black", bylevnames=c("White", "Black"), placement="topleft", add=TRUE, linecol=c("dark gray", "black"), secol="gray")
    plot(0:1, 0:1, type="n", ylab="", ylim=0:1, xlab="Sense of Linked Fate", main=paste(letters[5*(i-1)+4], ") Linked Fate\n", tab5b[i,4], sep=""), axes=FALSE)
    axis(1, at=seq(0, 1, 1/3), c("None", "Not much", "Some", "A lot"))
    axis(2, at=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), outcomeopts[[i]])
    axis(1, at=c(-99,99))
    axis(2, at=c(-99,99))
    axis(3, at=c(-99,99))
    axis(4, at=c(-99,99))
    abline(h=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), col="light gray", lty=3)
    plotwtdinteraction(linkedregs[[i]], "ml01", "black", bylevnames=c("White", "Black"), placement="topleft", add=TRUE, linecol=c("dark gray", "black"), secol="gray")
    plot(0:1, 0:1, type="n", ylab="", ylim=0:1, xlab="Group Salience", main=paste(letters[5*(i-1)+5], ") Group Salience\n", tab5b[i,5], sep=""), axes=FALSE)
    axis(1, at=seq(0, 1, 1/4), c("Lowest", "", "", "", "Highest"))
    axis(2, at=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), outcomeopts[[i]])
    axis(1, at=c(-99,99))
    axis(2, at=c(-99,99))
    axis(3, at=c(-99,99))
    axis(4, at=c(-99,99))
    abline(h=seq(0, 1, 1/(length(outcomeopts[[i]])-1)), col="light gray", lty=3)
    plotwtdinteraction(salienceregs[[i]], "salience", "black", bylevnames=c("White", "Black"), placement="topleft", add=TRUE, linecol=c("dark gray", "black"), secol="gray")
}
mtext("Perceptions by Race and", 3, outer=TRUE)
dev.off()



# Differences in Information-Seeking Preferences By Race

rcdifs <- as.data.frame(with(vpalt, t(sapply(baroutcomesreadchoice, function(x) unlist(wtd.t.test(nalevs(x)[black==FALSE], nalevs(x)[black==TRUE]))))), StringsAsFactors=FALSE)

rcdifstable <- data.frame(Whites=rcdifs[,"additional.Mean.x"], Blacks=rcdifs[,"additional.Mean.y"], Difference=rcdifs[,"additional.Difference"], sig=starmaker(as.numeric(as.character(rcdifs[,"coefficients.p.value"]))))

rownames(rcdifstable) <- rownames(rcdifs) <- c(paste("Excerpt", 1:8, c("Pro-Officer", "Pro-Taylor", "Pro-Officer", "Pro-Officer", "Pro-Officer", "Pro-Taylor", "Pro-Taylor", "Pro-Officer")), "All Pro-Officer", "All Pro-Taylor")

rcdifspaper <- rcdifstable[c("All Pro-Officer", "All Pro-Taylor"),]
write.csv(rcdifspaper, "Tables/TableH4-DesireToReadAdditional.csv")

rcdifs$padjust <- p.adjust(as.numeric(as.character(rcdifs$coefficients.p.value)), method="BY")

Nrcs <-  sapply(baroutcomesreadchoice, function(x) c(White=sum(!is.na(x) & vpalt$black==FALSE), Black=sum(!is.na(x) & vpalt$black==TRUE)))

rcomb <- data.frame(rcdifs[,c("additional.Mean.x", "additional.Mean.y", "additional.Difference", "additional.Std. Err", "coefficients.t.value", "coefficients.p.value", "padjust")], t(Nrcs))

rownames(rcomb) <- c(paste("Excerpt", 1:8, c("Pro-Officer", "Pro-Taylor", "Pro-Officer", "Pro-Officer", "Pro-Officer", "Pro-Taylor", "Pro-Taylor", "Pro-Officer")), "All Pro-Officer", "All Pro-Taylor")
colnames(rcomb) <- c("Whites Mean", "Blacks Mean", "Difference", "Standard Error", "t-statistic", "p-value", "adjusted p-value", "N of Whites", "N of Blacks")

write.csv(rcomb, "Tables/TableHG7-DesireToReadAllStatements.csv")

# Regressions for Difference in Information Seeking
readregs <- with(vpalt, c(lapply(baroutcomesreadchoice[1:8], function(y) glm(y~black+female+age01+ed01+pid01+libcon01, family="binomial")), lapply(baroutcomesreadchoice[9:10], function(y) lm(y~black+female+age01+ed01+pid01+libcon01))))

readregsNs <- sapply(readregs, function(x) length(predict(x)))
readregsR2s <- sapply(readregs, function(x) try(summary(x)$r.squared))
readregsPseudoR2s <- sapply(readregs, function(x) try(pR2(x)["McFadden"]))

readregsout <- rbind(rownames(rcomb)[1], ext(readregs[[1]]), N=readregsNs[1], R2=readregsR2s[1], PseudoR2=readregsPseudoR2s[1])
for(i in 2:length(w01regs)){
    regtemp <- rbind(rownames(rcomb)[i], ext(readregs[[i]]), N=readregsNs[i], R2=readregsR2s[i], PseudoR2=readregsPseudoR2s[i])
    readregsout <- cbind(readregsout, regtemp)
}

write.csv(readregsout, "Tables/TableHG8-DesireToReadAllStatementsRegressions.csv")



# Statement Weights Predicting Outcomes

appart <- nalevs(vpalt$approppart)

nopart <- lm(p01$appropriate~w01$pchiefweight+w01$s1weight+w01$s2weight+w01$s3weight+w01$s4weight+w01$s5weight+w01$s6weight, subset=!is.na(appart))
withpart <- lm(p01$appropriate~appart+w01$s3weight+w01$s4weight+w01$s5weight+w01$s6weight)
justpart <- lm(p01$appropriate~appart)


partwayanova <- anova(justpart, withpart)

chpart <- nalevs(vpalt$chargedpartway)

nopartcharge <- lm(p01$charged~w01$pchiefweight+w01$s1weight+w01$s2weight+w01$s3weight+w01$s4weight+w01$s5weight+w01$s6weight, subset=!is.na(chpart))
withpartcharge <- lm(p01$charged~chpart+w01$s3weight+w01$s4weight+w01$s5weight+w01$s6weight)

justpartcharge <- lm(p01$charged~chpart)

partchargeanova <- anova(justpartcharge, withpartcharge)



withpart2 <- with(w01, lm(p01$appropriate~appart*(s3weight+s4weight+s5weight+s6weight)))
withpartcharge2 <- with(w01, lm(p01$charged~appart*(s3weight+s4weight+s5weight+s6weight)))

postwt <- with(w01, (s3weight-s4weight-s5weight+s6weight))/2

withpart3 <- with(w01, lm(p01$appropriate~appart*postwt))
withpartcharge3 <- with(w01, lm(p01$charged~chpart*postwt))

pappk <- predict(withpart3, newdata=p01)
pchk <- predict(withpartcharge3, newdata=p01)

#lines(jitter(pappk[!is.na(appart)], amount=.05)~postwt[!is.na(appart)], type="p")
#lines(jitter(pchk[!is.na(chpart)], amount=.05)~postwt[!is.na(chpart)], type="p")

pdf("Figures/NoCutoff/Partial vs. Final Assessments By Relative Weights.pdf", width=16, height=8)
par(mfrow=c(1,2))
plot(c(-1,1.2), c(-.1,1.1), type="n", ylab="Final Appropriateness Assessment", xlab="Pro-Officer Weight - Pro-Victim Weight on subsequent Qs", axes=FALSE, main="Appropriateness as a function of partway assessments and relative\nweight given to pro-officer vs. pro-victim additional statments")
axis(1, c(-99,99))
axis(1, seq(-1,1,.2))
axis(2, c(-99,99))
axis(2, c(-99,99))
axis(2, seq(0, 1, .25), labels=gsub(" appropriate", "", levels(vpalt$approppart)))
axis(3, c(-99,99))
axis(4, c(-99,99))
abline(h=seq(0,1,.25), col="gray", lty=3)
plotwtdinteraction(withpart3, "postwt", "appart", add=TRUE, bylevnames=gsub(" appropriate", "", levels(vpalt$approppart)))
lines(jitter(pappk[!is.na(appart) & vpalt$black], amount=.02)~jitter(postwt[!is.na(appart) & vpalt$black], amount=.05), type="p", col="black", pch=21, bg="black", cex=.8)
lines(jitter(pappk[!is.na(appart) & !vpalt$black], amount=.02)~jitter(postwt[!is.na(appart) & !vpalt$black], amount=.05), type="p", col="black", pch=21, bg="gray", cex=.8)
plot(c(-1,1.2), c(-.1,1.1), type="n", ylab="Final Should Be Charged Assessment", xlab="Pro-Officer Weight - Pro-Victim Weight on subsequent Qs", axes=FALSE, main="Should be charged as a function of partway assessments and relative\nweight given to pro-officer vs. pro-victim additional statments")
axis(1, c(-99,99))
axis(1, seq(-1,1,.2))
axis(2, c(-99,99))
axis(2, c(-99,99))
axis(2, seq(0, 1, 1/3), labels=c("definitely NOT", "probably NOT", "probably", "definitely"))
axis(3, c(-99,99))
axis(4, c(-99,99))
abline(h=seq(0,1,1/3), col="gray", lty=3)
plotwtdinteraction(withpartcharge3, "postwt", "chpart", add=TRUE, bylevnames=c("definitely NOT", "probably NOT", "probably should", "definitely should"))
lines(jitter(pchk[!is.na(chpart) & vpalt$black], amount=.02)~jitter(postwt[!is.na(chpart) & vpalt$black], amount=.05), type="p", col="black", pch=21, bg="black", cex=.8)
lines(jitter(pchk[!is.na(chpart) & !vpalt$black], amount=.02)~jitter(postwt[!is.na(chpart) & !vpalt$black], amount=.05), type="p", col="black", pch=21, bg="gray", cex=.8)
dev.off()


pwnom <- c(names(coef(nopart))[1:4], "appart", "chpart", names(coef(nopart))[5:8])

partwaytable <- cbind(ext(nopart)[match(pwnom, rownames(ext(nopart))),], ext(withpart)[match(pwnom, rownames(ext(withpart))),], ext(nopartcharge)[match(pwnom, rownames(ext(nopartcharge))),], ext(withpartcharge)[match(pwnom, rownames(ext(withpartcharge))),])

partwaytablens <- c(length(predict(nopart)), "", "", "", length(predict(withpart)), "", "", "", length(predict(nopartcharge)), "", "", "", length(predict(justpartcharge)), "", "", "")

partwaytableR2s <- c(summary(nopart)$r.squared, "", "", "", summary(withpart)$r.squared, "", "", "", summary(nopartcharge)$r.squared, "", "", "", summary(justpartcharge)$r.squared, "", "", "")

partwaytableAnovaInfo <- c("", "", "", "", paste(rd(partwayanova[2,"F"], 1), "(", rd(partwayanova[2,"Df"], 0), ",", rd(partwayanova[2,"Res.Df"], 0), ")", sep=""), "", partwayanova[2,"Pr(>F)"], "", "", "", "", "",  paste(rd(partchargeanova[2,"F"], 1), "(", rd(partchargeanova[2,"Df"], 0), ",", rd(partchargeanova[2,"Res.Df"], 0), ")", sep=""), "", partchargeanova[2,"Pr(>F)"], "")

partwayexport <- rbind(partwaytable, n=partwaytablens, R2=partwaytableR2s, 'ANOVA F(df) p'=partwaytableAnovaInfo)

write.csv(cbind(partwayexport[,9:10], starmaker(as.numeric(as.character(partwayexport[,11]))), partwayexport[,13:14], starmaker(as.numeric(as.character(partwayexport[,15])))), "Tables/TableH5-PartWayAssessmentsChargedOnly.csv")

write.csv(partwayexport, "Tables/TableHG9-PartWayAssessmentsFull.csv")




table(appart, p01$appropriate)

appartchange <- p01$appropriate-appart

summary(lm(appartchange~vpalt$black))

#write.csv(coef(summary(lm(p01$appropriate~appart*vpalt$black))), "Tables/NoCutoff/PartwayChangesRegApprop.csv")



table(chpart, p01$charged)

chargechange <- p01$charged-chpart

summary(lm(chargechange~vpalt$black))

summary(lm(p01$charged~chpart*vpalt$black))

#write.csv(coef(summary(lm(p01$charged~chpart*vpalt$black))), "Tables/NoCutoff/PartwayChangesReCharged.csv")


black <- vpalt$black

findwtdinteraction(lm(chargechange~black), "black", acclevnames=c("White", "Black"))


appropchangeset <- table(cut(appartchange, c(-99, -.01, .01, 99), c("Less Appropriate", "Same", "More Appropriate")), vpalt$black)
chargedchangeset <- table(cut(chargechange, c(-99, -.01, .01, 99), c("Less Likely", "Same", "More Likely")), vpalt$black)
colnames(appropchangeset) <- colnames(chargedchangeset) <- c("White", "Black")

acs <- cbind(Overall=rowSums(appropchangeset), appropchangeset)
ccs <- cbind(Overall=rowSums(chargedchangeset), chargedchangeset)
    
appropchangetab <- rbind(t(t(acs)/colSums(acs)), n=colSums(acs))
chargechangetab <- rbind(t(t(ccs)/colSums(ccs)), n=colSums(ccs))

acouter <- rbind(appropchangetab, "", chargechangetab)
write.csv(acouter, "Tables/TableHK1-PartwayChanges.csv")




# Belief Updating With Controls

nopartcont <- with(vpalt, lm(p01$appropriate~w01$pchiefweight+w01$s1weight+w01$s2weight+w01$s3weight+w01$s4weight+w01$s5weight+w01$s6weight+black+female+age01+ed01+pid01+libcon01, subset=!is.na(appart)))
withpartcont <- with(vpalt, lm(p01$appropriate~appart+w01$s3weight+w01$s4weight+w01$s5weight+w01$s6weight+black+female+age01+ed01+pid01+libcon01))
justpartcont <- with(vpalt, lm(p01$appropriate~appart+black+female+age01+ed01+pid01+libcon01))

partwayanovacont <- anova(justpartcont, withpartcont)


nopartchargecont <- with(vpalt, lm(p01$charged~w01$pchiefweight+w01$s1weight+w01$s2weight+w01$s3weight+w01$s4weight+w01$s5weight+w01$s6weight+black+female+age01+ed01+pid01+libcon01, subset=!is.na(appart)))
withpartchargecont <- with(vpalt, lm(p01$charged~chpart+w01$s3weight+w01$s4weight+w01$s5weight+w01$s6weight+black+female+age01+ed01+pid01+libcon01))

justpartchargecont <- with(vpalt, lm(p01$charged~chpart+black+female+age01+ed01+pid01+libcon01))

partchargeanovacont <- anova(justpartchargecont, withpartchargecont)


pwnomcont <- c(names(coef(nopartcont))[1:4], "appart", "chpart", names(coef(nopartcont))[5:length(coef(nopartcont))])

partwaytablecont <- cbind(ext(nopartcont)[match(pwnomcont, rownames(ext(nopartcont))),], ext(withpartcont)[match(pwnomcont, rownames(ext(withpartcont))),], ext(nopartchargecont)[match(pwnomcont, rownames(ext(nopartchargecont))),], ext(withpartchargecont)[match(pwnomcont, rownames(ext(withpartchargecont))),])

rownames(partwaytablecont) <- pwnomcont

partwaytablecontns <- c(length(predict(nopartcont)), "", "", "", length(predict(withpartcont)), "", "", "", length(predict(nopartchargecont)), "", "", "", length(predict(justpartchargecont)), "", "", "")

partwaytablecontR2s <- c(summary(nopartcont)$r.squared, "", "", "", summary(withpartcont)$r.squared, "", "", "", summary(nopartchargecont)$r.squared, "", "", "", summary(justpartchargecont)$r.squared, "", "", "")

partwaytablecontAnovaInfo <- c("", "", "", "", paste(rd(partwayanovacont[2,"F"], 1), "(", rd(partwayanovacont[2,"Df"], 0), ",", rd(partwayanovacont[2,"Res.Df"], 0), ")", sep=""), "", partwayanovacont[2,"Pr(>F)"], "", "", "", "", "",  paste(rd(partchargeanovacont[2,"F"], 1), "(", rd(partchargeanovacont[2,"Df"], 0), ",", rd(partchargeanovacont[2,"Res.Df"], 0), ")", sep=""), "", partchargeanovacont[2,"Pr(>F)"], "")

partwaycontexport <- rbind(partwaytablecont, n=partwaytablecontns, R2=partwaytablecontR2s, 'ANOVA F(df) p'=partwaytablecontAnovaInfo)

write.csv(partwaycontexport, "Tables/TableHG10-PartWayAssessmentsFullwControls.csv")



# Tests of Reduction of Significance of Black Dummy

ocs <- data.frame(p01, w01)

#fd <- with(vpalt, data.frame(ocs, sexfac, ed01, age01, pid01, black, racres, differentialperceptions, salience))
#usable <- apply(fd, 1, function(x) sum(is.na(x))==0)

racebaseregs2 <- lapply(ocs[usable,], function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+pid01+libcon01+black)))
racepriorregs2 <- lapply(ocs[usable,], function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+pid01+libcon01+black+racres+differentialperceptions)))

vpalt$allprooff <- ocs$allprooff
vpalt$allprotaylor <- ocs$allprotaylor
vpalt$allprooffNoWalker <- ocs$allprooffNoWalker
dfset <- vpalt[usable,]

sensusable <- !is.na(rowSums(ocs)) & usable

outcomebaseregs <- with(vpalt[sensusable,], lapply(ocs[sensusable,], function(g) lm(g~sexfac+ed01+age01+pid01+libcon01+black)))
outcomefullregs <- with(vpalt[sensusable,], lapply(ocs[sensusable,], function(g) lm(g~sexfac+ed01+age01+pid01+libcon01+black+racres+differentialperceptions)))
racresfullregs <- with(vpalt[sensusable,], lm(racres~sexfac+ed01+age01+pid01+libcon01+black))
difpercfullregs <- with(vpalt[sensusable,], lm(differentialperceptions~sexfac+ed01+age01+pid01+libcon01+black))

cor(as.data.frame(sapply(outcomefullregs, resid)), data.frame(resid(racresfullregs), resid(difpercfullregs)))

write.csv(ext(difpercfullregs), "Tables/TableHG12APolicePred.csv")
write.csv(ext(racresfullregs), "Tables/TableHG12BRacResPred.csv")


bootn=1000
rbrboot1 <- lapply(1:bootn, function(x) sample(1:sum(usable), replace=TRUE))

raceeffects <- as.list(rep(NA, ncol(ocs)))
for(i in 1:ncol(ocs)){
    print(names(ocs)[i])
    newdf <- dfset[!is.na(ocs[,i]),]
    outcome <- newdf[,names(newdf)==names(ocs)[i]]
    basecoefs <- t(sapply(rbrboot1, function(g) coef(with(newdf[g,], lm(outcome[g]~sexfac+ed01+age01+pid01+libcon01+black)))))
    print("a")
    priorcoefs <- t(sapply(rbrboot1, function(g) coef(with(newdf[g,], lm(outcome[g]~sexfac+ed01+age01+pid01+libcon01+black+racres+differentialperceptions)))))
    print("b")
    racrescoefs <- t(sapply(rbrboot1, function(g) coef(with(newdf[g,], lm(racres~sexfac+ed01+age01+pid01+libcon01+black)))))
    print("c")
    difperccoefs <- t(sapply(rbrboot1, function(g) coef(with(newdf[g,], lm(differentialperceptions~sexfac+ed01+age01+pid01+libcon01+black)))))
    print("d")
    raceeffects[[i]] <- data.frame(total=basecoefs[,"blackTRUE"], direct=priorcoefs[,"blackTRUE"], indirectviaracres=racrescoefs[,"blackTRUE"]*priorcoefs[,"racres"], indirectviadifperc=difperccoefs[,"blackTRUE"]*priorcoefs[,"differentialperceptions"], totalindirect=(racrescoefs[,"blackTRUE"]*priorcoefs[,"racres"]+difperccoefs[,"blackTRUE"]*priorcoefs[,"differentialperceptions"]), racrescoef=priorcoefs[,"racres"], difperccoef=priorcoefs[,"differentialperceptions"])
}
names(raceeffects) <- names(ocs)

bootedcoefs <- lapply(raceeffects, function(x) t(apply(x, 2, function(g) c(mean=mean(g), se=sd(g), p=2*(1-max(table(g>0))/bootn)))))
propmediated <- sapply(raceeffects, function(x) mean(x$totalindirect/x$total))

write.csv(bootedcoefs, paste("Tables/TableHG11-FullBootstrapCoefficients.csv", sep=""))

collapsedboot <- t(sapply(bootedcoefs, function(x) apply(x, 1, function (g) paste(rd(g[1]), "(", rd(g[2]), ")", starmaker(g[3]), sep=""))))

paste(rd(ext(racresfullregs)["blackTRUE",][1]), "(", rd(ext(racresfullregs)["blackTRUE",][2]), ")", starmaker(ext(racresfullregs)["blackTRUE",][3]), sep="")

boottable <- cbind(collapsedboot[,c("total", "direct")], policebias=paste(rd(ext(difpercfullregs)["blackTRUE",][1]), "(", rd(ext(difpercfullregs)["blackTRUE",][2]), ")", starmaker(ext(difpercfullregs)["blackTRUE",][3]), sep=""), racres=paste(rd(ext(racresfullregs)["blackTRUE",][1]), "(", rd(ext(racresfullregs)["blackTRUE",][2]), ")", starmaker(ext(racresfullregs)["blackTRUE",][3]), sep=""), collapsedboot[,c("difperccoef", "racrescoef", "indirectviadifperc", "indirectviaracres", "totalindirect")])

write.csv(boottable, "Tables/TableH6-RaceEffectsMediations.csv")

##




shortocs <- ocs[usable,c("appropriate", "charged", "offattacked", "believeweapon", "racerole", "allprooff", "allprotaylor", "allprooffNoWalker")]

racepriorregsminbypart <- lapply(shortocs, function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+libcon01+(black+pid01)+(racres+differentialperceptions))))

racepriorregsminbypartinteract <- lapply(shortocs, function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+libcon01+(black*pid01)+(racres+differentialperceptions))))


racepriorregsbypart <- lapply(shortocs, function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+libcon01+(black+pid01)+(racres+differentialperceptions))))

racepriorregsbypartinteracts <- lapply(shortocs, function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+libcon01+(black*pid01)+(racres+differentialperceptions))))


racepriorregsbypartnointeractsbase <- lapply(shortocs, function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+libcon01+(black+pid01))))

racepriorregsbypartinteractsbase <- lapply(shortocs, function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+libcon01+(black*pid01))))

racepriorregsbypartinteracts2 <- lapply(shortocs, function(y) with(vpalt[usable,], lm(y~sexfac+ed01+age01+libcon01+(black*pid01)*(racres+differentialperceptions))))





interactanovas <- lapply(1:length(racepriorregsbypart), function(x) anova(racepriorregsbypart[[x]], racepriorregsbypartinteracts[[x]]))

interactanovasnocont <- lapply(1:length(racepriorregsbypartinteractsbase), function(x) anova(racepriorregsbypartnointeractsbase[[x]], racepriorregsbypartinteractsbase[[x]]))

iacadjp <- p.adjust(sapply(interactanovas, function(x) x[2,"Pr(>F)"]), "BY")
iacadjpnc <- p.adjust(sapply(interactanovasnocont, function(x) x[2,"Pr(>F)"]), "BY")

iacont <- paste(sapply(interactanovas, function(x) paste("F = ", rd(x[2,"F"], 1), "(", rd(x[2,"Df"], 0), "), p = ", rd(x[2,"Pr(>F)"]), sep="")), ", adj p = ", rd(iacadjp), sep="")

ianocont <- paste(sapply(interactanovasnocont, function(x) paste("F = ", rd(x[2,"F"], 1), "(", rd(x[2,"Df"], 0), "), p=", rd(x[2,"Pr(>F)"]), sep="")), ", adj p = ", rd(iacadjpnc), sep="")


shortocsnames <- c("Officer's Actions Appropriate", "Officer Should Be Charged", "Victim Attacked Officer", "Victim Had Weapon", "Race Role in Shooting", "Pro-Officer Statement Weights", "Pro-Victim Statement Weights", "Pro-Officer Statement Weights (Excluding Walker)")

son <- c("Officer's Actions\nAppropriate", "Officer Should\nBe Charged", "Victim Attacked\nOfficer", "Victim Had\nWeapon", "Race Role\nin Shooting", "Pro-Officer Statement\nWeights", "Pro-Victim Statement\nWeights", "Pro-Officer Statement\nWeights -Walker")



pdf("Figures/FigureHJ1-InfluenceOfRaceAndPartisanship.pdf", width=8, height=12)
par(mfrow=c(2,2), oma=c(0,3,6,0))
i <- 1
plotwtdinteraction(racepriorregsbypartinteractsbase[[i]], across="pid01", by="black", bylevnames=c("White", "Black"), acrosslevs=seq(0,1,1/6), acclevnames=c("Strong D", "", "", "Independent", "", "", "Strong R"), ylim=c(0,1), addat=TRUE, linecol=c("dark gray", "black"), secol=c("gray", "gray"), ylab="", main="", xlab="Partisanship", placement="topleft")
mtext("Influence of Race and Partisanship on Outcomes in Regression Models\n(F-Tests Measure Significance of Interactions)", outer=TRUE, line=2)
mtext(son[i], side=2, outer=FALSE, line=3.5)
mtext("No Controls for Racial Priors", outer=FALSE, line=2.5)
mtext(ianocont[i], outer=FALSE, line=1)
plotwtdinteraction(racepriorregsbypartinteracts[[i]], across="pid01", by="black", bylevnames=c("White", "Black"), acrosslevs=seq(0,1,1/6), acclevnames=c("Strong D", "", "", "Independent", "", "", "Strong R"), ylim=c(0,1), addat=TRUE, linecol=c("dark gray", "black"), secol=c("gray", "gray"), ylab="", main="", xlab="Partisanship", placement="topleft")
mtext("Controlling for Racial Priors", outer=FALSE, line=2.5)
mtext(iacont[i], outer=FALSE, line=1)
for(i in 2:length(racepriorregsbypartinteracts)){
mtext("Influence of Race and Partisanship on Outcomes in Regression Models\n(F-Tests Measure Significance of Interactions)", outer=TRUE, line=2)
plotwtdinteraction(racepriorregsbypartinteractsbase[[i]], across="pid01", by="black", bylevnames=c("White", "Black"), acrosslevs=seq(0,1,1/6), acclevnames=c("Strong D", "", "", "Independent", "", "", "Strong R"), ylim=c(0,1), addat=TRUE, linecol=c("dark gray", "black"), secol=c("gray", "gray"), ylab="", main="", xlab="Partisanship", placement="topleft")
mtext(son[i], side=2, outer=FALSE, line=3.5)
mtext("No Controls for Racial Priors", outer=FALSE, line=2.5)
mtext(ianocont[i], outer=FALSE, line=1)
plotwtdinteraction(racepriorregsbypartinteracts[[i]], across="pid01", by="black", bylevnames=c("White", "Black"), acrosslevs=seq(0,1,1/6), acclevnames=c("Strong D", "", "", "Independent", "", "", "Strong R"), ylim=c(0,1), addat=TRUE, linecol=c("dark gray", "black"), secol=c("gray", "gray"), ylab="", main="", xlab="Partisanship", placement="topleft")
mtext("Controlling for Racial Priors", outer=FALSE, line=2.5)
mtext(iacont[i], outer=FALSE, line=1)
}
dev.off()





pdf("Figures/FigureHJ2-InfluenceOfPoliceExpectationsByRaceAndPartisanship(Omitted).pdf", width=8, height=12)
par(mfrow=c(2,2), oma=c(0,3,6,0))
for(i in 1:length(racepriorregsbypartinteracts)){
plotwtdinteraction(racepriorregsbypartinteracts[[i]], across="differentialperceptions", by="black", at="pid01", atlevs=c(0,1), acrosslevs=seq(.3, 1, .1), bylevnames=c("White", "Black"), atlevnames=c("Strong D", "Strong R"), ylim=c(0,1), addat=TRUE, linecol=c("Blue", "Red"), secol=c("light blue", "pink"), ylab=shortocsnames[i], main="Police Fairness", xlab="Police Bias")#paste("Interaction Plot of", shortocsnames[i], "Across\nPerceptions of Police Fairness By Race and Partisanship")
plotwtdinteraction(racepriorregsbypartinteracts[[i]], across="racres", by="black", at="pid01", atlevs=c(0,1), acrosslevs=seq(0, 1, .1), bylevnames=c("White", "Black"), atlevnames=c("Strong D", "Strong R"), ylim=c(0,1), addat=TRUE, linecol=c("Blue", "Red"), secol=c("light blue", "pink"), ylab=shortocsnames[i], main="Racial Resentment", xlab="Racial Resentment")#paste("Interaction Plot of", shortocsnames[i], "Across\nRacial Resentment By Race and Partisanship")
mtext("Interaction Plots of Outcomes by Race and Partisanship", outer=TRUE, line=2.5)
}
dev.off()

pdf("Figures/FigureHJ3-InfluenceOfPoliceExpectationsInteractedWithRaceAndPartisanship(Omitted).pdf", width=8, height=12)
par(mfrow=c(2,2), oma=c(0,3,6,0))
for(i in 1:length(racepriorregsbypartinteracts2)){
plotwtdinteraction(racepriorregsbypartinteracts2[[i]], across="differentialperceptions", by="black", at="pid01", atlevs=c(0,1), acrosslevs=seq(.3, 1, .1), bylevnames=c("White", "Black"), atlevnames=c("Strong D", "Strong R"), ylim=c(0,1), addat=TRUE, linecol=c("Blue", "Red"), secol=c("light blue", "pink"), ylab=shortocsnames[i], main="Police Fairness", xlab="Police Bias")#paste("Interaction Plot of", shortocsnames[i], "Across\nPerceptions of Police Fairness By Race and Partisanship")
plotwtdinteraction(racepriorregsbypartinteracts2[[i]], across="racres", by="black", at="pid01", atlevs=c(0,1), acrosslevs=seq(0, 1, .1), bylevnames=c("White", "Black"), atlevnames=c("Strong D", "Strong R"), ylim=c(0,1), addat=TRUE, linecol=c("Blue", "Red"), secol=c("light blue", "pink"), ylab=shortocsnames[i], main="Racial Resentment", xlab="Racial Resentment")#paste("Interaction Plot of", shortocsnames[i], "Across\nRacial Resentment By Race and Partisanship")
mtext("Interaction Plots of Outcomes by Race and Partisanship", outer=TRUE, line=2.5)
}
dev.off()


interactanovas <- lapply(1:length(racepriorregsbypart), function(x) anova(racepriorregsbypart[[x]], racepriorregsbypartinteracts[[x]]))

interactanovas2 <- lapply(1:length(racepriorregsbypart), function(x) anova(racepriorregsminbypart[[x]], racepriorregsminbypartinteract[[x]]))



eachpredoutcome <- with(vpalt, lapply(p01, function(y) t(sapply(seppredRRPPs, function(x) coef(summary(lm(y~x+sexfac+ed01+age01+libcon01+black+pid01)))[2,]))))

write.csv(as.data.frame(eachpredoutcome), "Tables/TableHY1-SeparatePredictorsFromScales.csv")





library(vioplot)
pdf("Figures/FigureHG2-BlackWhiteDifferences.pdf", width=12, height=8)
par(mfrow=c(1,2))
with(vpalt, vioplot(racres~whiteblack, ylim=0:1, xlab="Race", ylab="Racial Resentment Score", main="Racial Differences in Racial Resentment"))
with(vpalt, vioplot(differentialperceptions~whiteblack, ylim=0:1, xlab="Race", ylab="Perceptions of Police Bias", main="Racial Differences in Perceptions of Police Bias"))
dev.off()


with(vpalt, vioplot(salience[black]~condition[black], ylim=0:1, xlab="Condition"))
with(vpalt, vioplot(salience[!black]~condition[!black], ylim=0:1, xlab="Condition"))

